#include<stdio.h>
#include<math.h>
#include<stdlib.h>


#define GAP 100
#define PART 20.0
#define INIT_WT_RANGE 0.5
#define L_RATEW_LIMIT 0.2
#define INIT_MULT -5.0

/***VARIABLE DECLARATIONS****/


int *gia1HidNodesNo;
int giDim;
int giClassNo;
int giTotalNodes;
int giTrainDataCount;
int giTestDataCount;
int *gia1TrainTarget;
int *gia1TestTarget;
int giFinItrCount;
int giCurItrCount;
int giMSTE;
int giMSTR;

double gda1MaxMin[2];
double ***gda3w;
double *gda1mult;
double *gda1new_mult;
double *gda1InputVec;
double *gda1TargetVec;
double **gda2OutVec;
double **gda2DerOut;
double **gda2TrainData;
double **gda2TestData;
double **del;
double gdLRateW;
double gdCurMSETrainSet;
double gdCurMSETestSet;
double gdPrevMSETrainSet;
double gdPrevMSETestSet;
double gdMuRate=0.5;
double *gda1DerMult;

/*****FUNCTION DECLARATIONS***/


void input(void);

	void take_architecture(void);
	void allocate_memory(void);
	void put_data_into_array(void);
	void normalize(void);

		void max_find(void);
		void alter_data(void);

	void set_parameters(void);

void initialise(void);

	void initialise_network_weights(void);

void train(void);

	void output_evaluate(void);

		void out_zero(void);
		void display_performance(void);

	void balance_wts(void);
	void performance_analysis(void);

		void find_MSE_train_set(void);
		void find_MSE_test_set(void);

			double err(void);

		void find_MSTR(void);
		void find_MSTE(void);


void the_end(void);


main()
{
input();
initialise();
train();
the_end();
return;
}

void input(void)
{
take_architecture();
allocate_memory();
put_data_into_array();
normalize();
set_parameters();
return;
}

void take_architecture(void)
{

int HidNo,i;
char c;

printf("\nenter in the number of features:- ");
scanf("%d",&giDim);fflush(stdin);
printf("\nenter in the no. of classes:- ");
scanf("%d",&giClassNo);fflush(stdin);
printf("\nenter in the no. of hidden layers:- ");
scanf("%d",&HidNo);fflush(stdin);

if((gia1HidNodesNo=(int*)malloc((HidNo+2)*sizeof(int)))==NULL)
printf("Error Allocating Memory for gia1HidNodesNo\n");
gia1HidNodesNo[0]=HidNo;

for(i=1;i<=gia1HidNodesNo[0];i++)
 {
	printf("\nenter in the no. of nurons in layer no.%d:- ",i);
	scanf("%d",&gia1HidNodesNo[i]);fflush(stdin);
	giTotalNodes+=gia1HidNodesNo[i];
 }

gia1HidNodesNo[gia1HidNodesNo[0]+1]=giClassNo;

if((gda1mult = (double *)malloc((giDim+1)*sizeof(double)))==NULL)
printf("Error Allocating Memory for gda1multo\n");
if((gda1new_mult = (double *)malloc((giDim+1)*sizeof(double)))==NULL)
printf("Error Allocating Memory for gda1new_mult\n");
if((gda1DerMult = (double *)malloc((giDim+1)*sizeof(double)))==NULL)
printf("Error Allocating Memory for gda1DerMult\n");
for(i=1;i<=giDim;i++)
  {
  gda1mult[i] = INIT_MULT;
  gda1new_mult[i] = INIT_MULT;
  }
/*printf("\n are you satisfied(y/n):- ");
//scanf("%c",&c);fflush(stdin);
//if(c=='n') take_architecture();
*/
return;
}

void allocate_memory(void)
{

int i,j;

if((gda3w=(double***)malloc((gia1HidNodesNo[0]+1)*sizeof(double**)))==NULL)
printf("Error Allocating Memory for gda3w\n");
if((gda3w[0]=(double**)malloc((gia1HidNodesNo[1]+1)*sizeof(double*)))==NULL)
printf("Error Allocating Memory for gda3w[0]\n");

for(i=0;i<=gia1HidNodesNo[1];i++)
   {
	if((gda3w[0][i]=(double*)malloc((giDim+1)*sizeof(double)))==NULL)
		printf("Error Allocating Memory for gda3w[0][i]\n");
   }

for(i=1;i<=gia1HidNodesNo[0];i++)
   {

	if((gda3w[i]=(double**)malloc((gia1HidNodesNo[i+1]+1)*sizeof(double*)))==NULL)
		printf("Error Allocating Memory for gda3w[i]\n");
	for(j=0;j<=gia1HidNodesNo[i+1];j++)
	 {
		if((gda3w[i][j]=(double*)malloc((gia1HidNodesNo[i]+1)*sizeof(double)))==NULL)
			printf("Error Allocating Memory for gda3w[i][j]\n");
	 }

   }

if((gda1InputVec=(double*)malloc((giDim+1)*sizeof(double)))==NULL)
			printf("Error Allocating Memory for gda1InputVec\n");
if((gda1TargetVec=(double*)malloc((giClassNo+1)*sizeof(double)))==NULL)
			printf("Error Allocating Memory for gda1TargetVec\n");
if((gda2OutVec=(double**)malloc((gia1HidNodesNo[0]+2)*sizeof(double*)))==NULL)
			printf("Error Allocating Memory for gda2OutVec\n");

for(i=0;i<=gia1HidNodesNo[0]+1;i++){
	if((gda2OutVec[i]=(double*)malloc((gia1HidNodesNo[i]+1)*sizeof(double)))==NULL)
		printf("Error Allocating Memory for gda2OutVec[i]\n");
}
if((gda2DerOut=(double**)malloc((gia1HidNodesNo[0]+2)*sizeof(double*)))==NULL)
		printf("Error Allocating Memory for gda2DerOut\n");

for(i=0;i<=gia1HidNodesNo[0]+1;i++){
	if((gda2DerOut[i]=(double*)malloc((gia1HidNodesNo[i]+1)*sizeof(double)))==NULL)
		printf("Error Allocating Memory for gda2DerOut[i]\n");
}
if((del=(double**)malloc((gia1HidNodesNo[0]+2)*sizeof(double*)))==NULL)
		printf("Error Allocating Memory for del\n");
for(i=0;i<=gia1HidNodesNo[0]+1;i++){
	if((del[i]=(double*)malloc((gia1HidNodesNo[i]+1)*sizeof(double)))==NULL)
		printf("Error Allocating Memory for del[i]\n");
}

return;
}


void put_data_into_array(void)
{
int i,j,pond1;
char file1[30];
FILE *filer;
double pond;

printf("\ngive in the name of your training data file:- ");
scanf("%s",file1);fflush(stdin);
filer=fopen(file1,"r");
giTrainDataCount=0;
giTestDataCount=0;
while(!feof(filer))
  {
	for(i=1;i<=giDim;i++)
	fscanf(filer,"%lf",&pond);
	fscanf(filer,"%d",&pond1);
	fscanf(filer,"\n");
	giTrainDataCount++;
  }

rewind(filer);

if((gda2TrainData=(double**)malloc(giTrainDataCount*sizeof(double*)))==NULL)
		printf("Error Allocating Memory for gda2TrainData\n");
for(i=0;i<giTrainDataCount;i++){
	if((gda2TrainData[i]=(double*)malloc((giDim+1)*sizeof(double)))==NULL)
				printf("Error Allocating Memory for gda2TrainData[i]\n");
}
if((gia1TrainTarget=(int*)malloc(giTrainDataCount*sizeof(int)))==NULL)
				printf("Error Allocating Memory for gia1TrainTarget\n");

for(j=0;j<giTrainDataCount;j++)
  {
	for(i=1;i<=giDim;i++)
		fscanf(filer,"%lf",&gda2TrainData[j][i]);
	fscanf(filer,"%d",&gia1TrainTarget[j]);
	fscanf(filer,"\n");

/* Patch by NRP to so that labels numbered like 1, 2, 3 can
be used **/
gia1TrainTarget[j]--;
  }
fclose(filer);

printf("\ngive in the name of your testing data file:- ");
scanf("%s",file1);fflush(stdin);
filer=fopen(file1,"r");
while(!feof(filer))
  {
	for(i=1;i<=giDim;i++)
		fscanf(filer,"%lf",&pond);
	fscanf(filer,"%d",&pond1);
	fscanf(filer,"\n");
	giTestDataCount++;
  }
rewind(filer);


printf("no of test data points %d \n", giTestDataCount);
printf("no of training data %d \n",giTrainDataCount);

if((gda2TestData=(double**)malloc(giTestDataCount*sizeof(double*)))==NULL)
				printf("Error Allocating Memory for gda2TestData\n");

for(i=0;i<giTestDataCount;i++){
	if((gda2TestData[i]=(double*)malloc((giDim+1)*sizeof(double)))==NULL)
				printf("Error Allocating Memory for gda2TestData[i]\n");
}

if((gia1TestTarget=(int*)malloc(giTestDataCount*sizeof(int)))==NULL)
				printf("Error Allocating Memory for gia1TestTarget\n");
for(j=0;j<giTestDataCount;j++)
  {
	for(i=1;i<=giDim;i++)
		fscanf(filer,"%lf",&gda2TestData[j][i]);
	fscanf(filer,"%d",&gia1TestTarget[j]);
	fscanf(filer,"\n");

/* Patch by NRP to so that labels numbered like 1, 2, 3 can
be used **/
gia1TestTarget[j]--;

  }
fclose(filer);

return;
}

void normalize(void)
{
max_find();
alter_data();
return;
}

void max_find(void)
{

int i,k;

gda1MaxMin[0]=fabs(gda2TrainData[0][1]);
gda1MaxMin[1]=gda1MaxMin[0];

for(i=2;i<=giDim;i++)
 {
	if(fabs(gda2TrainData[0][i]) > gda1MaxMin[1]) gda1MaxMin[1] = fabs(gda2TrainData[0][i]);
	if(fabs(gda2TrainData[0][i]) < gda1MaxMin[0]) gda1MaxMin[0] = fabs(gda2TrainData[0][i]);
 }

for(k=1;k<giTrainDataCount;k++)
 {
	for(i=1;i<=giDim;i++)
	  {
		if(fabs(gda2TrainData[k][i]) > gda1MaxMin[1]) gda1MaxMin[1] = fabs(gda2TrainData[k][i]);
		if(fabs(gda2TrainData[k][i]) < gda1MaxMin[0]) gda1MaxMin[0] = fabs(gda2TrainData[k][i]);
	  }
 }

for(k=0;k<giTestDataCount;k++)
 {
	for(i=1;i<=giDim;i++)
	  {
		if(fabs(gda2TestData[k][i]) > gda1MaxMin[1]) gda1MaxMin[1] = fabs(gda2TestData[k][i]);
		if(fabs(gda2TestData[k][i]) < gda1MaxMin[0]) gda1MaxMin[0] = fabs(gda2TestData[k][i]);
	  }
 }

return;
}



void alter_data(void)
{

int i,k;

for(k=0;k<giTrainDataCount;k++)
 {
	for(i=1;i<=giDim;i++)
            gda2TrainData[k][i]=gda2TrainData[k][i]/gda1MaxMin[1];
 }
for(k=0;k<giTestDataCount;k++)
 {
	for(i=1;i<=giDim;i++)
            gda2TestData[k][i]=gda2TestData[k][i]/gda1MaxMin[1];
 }

return;
}


void set_parameters(void)
{

printf("\n Enter ni the number of itterations:-");
scanf("%d",&giFinItrCount);fflush(stdin);
printf("\ngive in the learning rate for the weights:- ");
scanf("%lf",&gdLRateW);fflush(stdin);
printf("\n Enter the multiplier rate:- ");
scanf("%lf",&gdMuRate);
/*multi=lm;ni=n;*/

return;
}

void initialise(void)
{
initialise_network_weights();
return;
}

void initialise_network_weights(void)
{
/*unsigned*/
long seed;
int i,j,k;

printf("\n Enter in a seed for initialization:-");
scanf("%ld",&seed);fflush(stdin);
srand48(seed);
for(i=1;i<=gia1HidNodesNo[1];i++)
 {
	for(j=1;j<=giDim;j++)
	  {
            gda3w[0][i][j]=INIT_WT_RANGE*(-1.0+2*drand48());
	  }
  }

for(i=1;i<=gia1HidNodesNo[0];i++)
  {
	for(j=1;j<=gia1HidNodesNo[i+1];j++)

	  {
		for(k=1;k<=gia1HidNodesNo[i];k++)

                   gda3w[i][j][k]=INIT_WT_RANGE*(-1.0+2*drand48());
	  }
   }

return;
}


void train(void)
{
double l=gdLRateW;
int i,k,p1=0,n;
char cc='y';

giCurItrCount=0;

gdCurMSETrainSet=0;
gdCurMSETestSet=0;


performance_analysis();
display_performance();
gdPrevMSETrainSet=gdCurMSETrainSet;
gdPrevMSETestSet=gdCurMSETestSet;

while(cc=='y')
  {
	while(giCurItrCount < giFinItrCount )
	  {


		if(gdLRateW > L_RATEW_LIMIT)
		  {
			if(gdCurMSETrainSet > gdPrevMSETrainSet)
			  {
				p1++;

				if(p1 > 5)
				  {
					p1=0;
					gdLRateW=gdLRateW-l/PART;
				  }
			  }

			else p1=0;
		  }

		if(l < L_RATEW_LIMIT) l=L_RATEW_LIMIT;

		for(k=0;k<giTrainDataCount;k++)
		  {
			for(i=1;i<=giDim;i++)
				gda1InputVec[i]=gda2TrainData[k][i];

			for(i=1;i<=giClassNo;i++)
			  {
				if(i==gia1TrainTarget[k]+1) gda1TargetVec[i]=1.0;
				else gda1TargetVec[i]=0.0;
			  }

			output_evaluate();
			balance_wts();
		  }

		giCurItrCount++;
		performance_analysis();

		if(giCurItrCount % GAP == 0)
		  {
			display_performance();
                        printf("\n");
                        for(i=1;i<=giDim;i++)
			  {
                            printf("%lf ",gda1mult[i]);
			  }
			gdPrevMSETrainSet=gdCurMSETrainSet;
			gdPrevMSETestSet=gdCurMSETestSet;
		  }
	    }
	printf("\n How many more iterations:- ");
	scanf("%d",&n);
	if(n==0)
	  {
		cc='n';goto pip;
	  }
	giFinItrCount+=n;
	pip:;
  }
return;
}

void output_evaluate(void)
{
   int i,j,k;
   double att;
   for(i=1;i<=gia1HidNodesNo[0]+1;i++)
   {
    for(j=1;j<=gia1HidNodesNo[i];j++)
    {
    gda2OutVec[i][j]=0.0;
    }
   }
   for(i=1;i<=giDim;i++)
     {
       att = 1.0/(1.0+exp(-1.0*gda1mult[i]));
       gda1DerMult[i] = att*(1.0-att); 
     }
   for(i=1;i<=gia1HidNodesNo[1];i++)
   {
    for(j=1;j<=giDim;j++)
    {
    att = 1.0/(1.0+exp(-1.0*gda1mult[j])); 
    gda2OutVec[1][i]=gda2OutVec[1][i]+gda3w[0][i][j]*gda1InputVec[j]*att;
    }
    gda2OutVec[1][i]=1.0/(1.0+exp(-1.0*gda2OutVec[1][i]));
    gda2DerOut[1][i]=gda2OutVec[1][i]*(1.0-gda2OutVec[1][i]);
   }
   for(i=1;i<=gia1HidNodesNo[0];i++)
   {
    for(j=1;j<=gia1HidNodesNo[i+1];j++)
    {
     for(k=1;k<=gia1HidNodesNo[i];k++)
     {
     gda2OutVec[i+1][j]=gda2OutVec[i+1][j]+gda2OutVec[i][k]*gda3w[i][j][k];
     }
     gda2OutVec[i+1][j]=1.0/(1.0+exp(-1.0*gda2OutVec[i+1][j]));
     gda2DerOut[i+1][j]=gda2OutVec[i+1][j]*(1.0-gda2OutVec[i+1][j]);
    }}
}

void balance_wts(void)
{
int i,j,k;
double sum;
double att;

/***DELS AND DEL1 ARE BEING INITIALIZED*****/

for(i=1;i<=gia1HidNodesNo[0]+1;i++)
 {
	for(j=0;j<=gia1HidNodesNo[i];j++)
	  {
		del[i][j]=0.0;
	  }
 }

for(i=1;i<=giClassNo;i++)
	del[gia1HidNodesNo[0]+1][i]=gda2DerOut[gia1HidNodesNo[0]+1][i]*(gda1TargetVec[i]-gda2OutVec[gia1HidNodesNo[0]+1][i]);
     /***calculating deltas for all hidden layers *****/

for(i=gia1HidNodesNo[0];i>=1;i--)
 {
       for(k=1;k<=gia1HidNodesNo[i];k++)
	 {
		for(j=1;j<=gia1HidNodesNo[i+1];j++)
			del[i][k]=del[i][k]+del[i+1][j]*gda3w[i][j][k];		               del[i][k]=gda2DerOut[i][k]*del[i][k];
	 }

 }

for(i=1;i<=giDim;i++)
  {
    del[0][i]=0;
    for(j=0;j<=gia1HidNodesNo[1];j++)
      {
       del[0][i] =del[0][i] + del[1][j]*gda3w[0][j][i];
      }
    del[0][i] = del[0][i]*gda1DerMult[i];
  }
 
     /***balancing weights from output layer to last but one****/

for(i=1;i<=gia1HidNodesNo[0];i++)
 {
	for(j=1;j<=gia1HidNodesNo[i+1];j++)
	  {
		for(k=1;k<=gia1HidNodesNo[i];k++)
			gda3w[i][j][k]=gdLRateW*del[i+1][j]*gda2OutVec[i][k]+gda3w[i][j][k];
	  }
 }


/****balancing for the weights from input to 1 hidden layer*****/


for(j=1;j<=giDim;j++)
{
 att = 1.0/(1.0+exp(-1.0*gda1mult[j]));
for(i=1;i<=gia1HidNodesNo[1];i++)
 {
    gda3w[0][i][j]=gdLRateW*del[1][i]*gda1InputVec[j]*att+gda3w[0][i][j];
          }
 }



/***balancing the multipliers *****/

/*** for(i=1;i<giDim;i++) ***/


for(i=1;i<=giDim;i++)
      gda1mult[i] = gda1mult[i] + gdMuRate*gda1InputVec[i]*del[0][i];

return;
}

void performance_analysis(void)
{
find_MSE_train_set();
if(giCurItrCount % GAP == 0)
  {
	find_MSE_test_set();
	find_MSTR();
	find_MSTE();
  }
return;
}

void find_MSE_train_set(void)
{
int i,k;
double error=0;


for(k=0;k<giTrainDataCount;k++)
  {
	for(i=1;i<=giDim;i++)
		gda1InputVec[i]=gda2TrainData[k][i];

	for(i=1;i<=giClassNo;i++)
	  {
		if(i==gia1TrainTarget[k]+1)
			gda1TargetVec[i]=1.0;
		else
			gda1TargetVec[i]=0.0;
	  }
	output_evaluate();
	error+=err();
  }

error=error/((double)(giTrainDataCount*giClassNo));
gdCurMSETrainSet=error;

return;
}

void find_MSE_test_set(void)
{
int i,k;
double error=0;

for(k=0;k<giTestDataCount;k++)
  {
	for(i=1;i<=giDim;i++)
		gda1InputVec[i]=gda2TestData[k][i];

	for(i=1;i<=giClassNo;i++)
	  {
		if(i==gia1TestTarget[k]+1) gda1TargetVec[i]=1.0;
		else gda1TargetVec[i]=0.0;
	  }
	output_evaluate();
	error+=err();
  }

error=error/((double)(giTestDataCount*giClassNo));
gdCurMSETestSet=error;
return;
}

double err(void)
{
double error=0;
int i;

for(i=1;i<=giClassNo;i++)
 {
	error+=(gda1TargetVec[i]-gda2OutVec[gia1HidNodesNo[0]+1][i])*(gda1TargetVec[i]-gda2OutVec[gia1HidNodesNo[0]+1][i]);
 }

return(error);
}


void find_MSTR(void)
{
int i,j,k;
int cnt=0;
double val;

for(k=0;k<giTrainDataCount;k++)
  {
	for(i=1;i<=giDim;i++)
		gda1InputVec[i]=gda2TrainData[k][i];

	for(i=1;i<=giClassNo;i++)
	  {
		if(i==gia1TrainTarget[k]+1) gda1TargetVec[i]=1.0;
		else gda1TargetVec[i]=0.0;
	  }
	output_evaluate();

	j=1;
	val=gda2OutVec[gia1HidNodesNo[0]+1][1];
	for(i=2;i<=giClassNo;i++)
	  {
		if(val<gda2OutVec[gia1HidNodesNo[0]+1][i])
		  {
			val=gda2OutVec[gia1HidNodesNo[0]+1][i];
			j=i;
		  }
	  }
	if(j!=gia1TrainTarget[k]+1) cnt++;
  }

giMSTR=cnt;
return;
}

void find_MSTE(void)
{
int i,j,k;
int cnt=0;
double val;

for(k=0;k<giTestDataCount;k++)
  {
	for(i=1;i<=giDim;i++)
		gda1InputVec[i]=gda2TestData[k][i];

	for(i=1;i<=giClassNo;i++)
	  {
		if(i==gia1TestTarget[k]+1) gda1TargetVec[i]=1.0;
		else gda1TargetVec[i]=0.0;
	  }
	output_evaluate();

	j=1;
	val=gda2OutVec[gia1HidNodesNo[0]+1][1];
	for(i=2;i<=giClassNo;i++)
	  {
		if(val<gda2OutVec[gia1HidNodesNo[0]+1][i])
		  {
			val=gda2OutVec[gia1HidNodesNo[0]+1][i];
			j=i;
		  }
	  }
	if(j!=gia1TestTarget[k]+1) cnt++;
  }

giMSTE=cnt;

return;
}

void display_performance(void)
{

printf("\nittrn= %d : MSEtr=%lf, MSEte=%lf, MSTR=%d, MSTE=%d",giCurItrCount,gdCurMSETrainSet,gdCurMSETestSet,giMSTR,giMSTE);
return;
}

void the_end(void)
{

return;
}
















